Fixing unit conversion for C-grid velocity interpolation #2455

erikvansebille · 2026-01-08T15:59:46Z

This PR reverts back the change in c311fba, to avoid having to set fieldset.V.units = GeographicPolar()

Note that this aligns with the implementation of unit conversion in v3, where u and v are converted by the same meshJac

Parcels/parcels/field.py

Lines 1615 to 1635 in 1ce5908

    
           rad = np.pi / 180.0 
        
           deg2m = 1852 * 60.0 
        
           if applyConversion: 
        
               meshJac = (deg2m * deg2m * math.cos(rad * y)) if grid.mesh == "spherical" else 1 
        
           else: 
        
               meshJac = deg2m if grid.mesh == "spherical" else 1 
        
           jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) * meshJac 
        
           u = ( 
        
               (-(1 - eta) * U - (1 - xsi) * V) * px[0] 
        
               + ((1 - eta) * U - xsi * V) * px[1] 
        
               + (eta * U + xsi * V) * px[2] 
        
               + (-eta * U + (1 - xsi) * V) * px[3] 
        
           ) / jac 
        
           v = ( 
        
               (-(1 - eta) * U - (1 - xsi) * V) * py[0] 
        
               + ((1 - eta) * U - xsi * V) * py[1] 
        
               + (eta * U + xsi * V) * py[2] 
        
               + (-eta * U + (1 - xsi) * V) * py[3] 
        
           ) / jac

This has been the case since the very first implementation of C-grid interpolation in 19ef089

The question, however, is whether this is indeed correct to do. I never noticed it, but it may feel a bit strange to apply the same unit converter to lon and lat.

From discussions today, I now realise that our 'gold standard test', the zonal flow in nemo_curvilinear, in not a good test for this case, because v==np.close(0.) in that case, so it doesn't matter if it's decided by 1852*60 or by 1852*60*cos(lat)

I will go back to Philippe's original paper to see why we originally implemented the same conversion for u and v in Parcels v2.0. Hopefully that clarifies this....

To test that the spherical mesh conversion leads to the expected meridional velocity

erikvansebille · 2026-01-09T08:19:05Z

OK, I dug a bit deeper and now think that the original interpolation in v3 (and v4) is correct. See the new unit test I created in 5436136.

Parcels/tests/test_advection.py

Lines 459 to 481 in 5436136

    
           def test_nemo_curvilinear_with_vvel(): 
        
               """Testing that a constant meridional velocity field in a NEMO curvilinear grid 
        
               results in correct particle movement. 
        
               Only works in Southern hemisphere (away from the tripolar rotated grid). 
        
               """ 
        
               data_folder = parcels.download_example_dataset("NemoCurvilinear_data") 
        
               ds_fields = xr.open_mfdataset( 
        
                   data_folder.glob("*.nc4"), 
        
                   data_vars="minimal", 
        
                   coords="minimal", 
        
                   compat="override", 
        
               ) 
        
               ds_fields = ds_fields.isel(time=0, z_a=0, z=0, drop=True) 
        
               vvel = 0.5  # m/s 
        
               ds_fields["V"].data[:] = vvel  # set a constant meridional velocity 
        
               fieldset = parcels.FieldSet.from_nemo(ds_fields) 
        
               lonp, latp = 30, -70 
        
               pset = parcels.ParticleSet(fieldset, lon=lonp, lat=latp) 
        
               pset.execute(AdvectionEE, runtime=np.timedelta64(10, "D"), dt=np.timedelta64(1, "D")) 
        
               v = (pset.lat - latp) / pset.time * 1852 * 60 
        
               np.testing.assert_allclose(v, vvel, atol=1e-5)

This sets a meridional vvel=0.5 on all grid points of the fieldset.V of NemoCurvilinear_data, and then tests whether the average velocity of particles advected indeed is vvel (accounting for the fact that one degree of latitude is 1852*60 meters).

The test passes with the current v4 (and also v3; I checked manually), but not with the change below in /src/parcels/_core/field.py

@@ -312,8 +312,8 @@ class VectorField:
 
         if applyConversion:
             if self.U.grid._mesh == "spherical":
-                meshJac = 1852 * 60.0 * np.cos(np.deg2rad(y))
-                u = u / meshJac
+                meshJac = 1852 * 60.0
+                u = u / (meshJac * np.cos(np.deg2rad(y)))
                 v = v / meshJac
             if "3D" in self.vector_type:
                 w = self.W.units.to_target(w, z, y, x)

So this means that u and v need to be converted with latitude-dependence. I find it counter-intuitive (anyone who can explain me why; please be my guest!), but the good news is that this does not seem to be a bug ☺️

MeikeBos · 2026-01-09T11:06:25Z

I think the same unit converter meshJac for both u and v arises from the fact that we use the fluxes to calculate the velocities for the c-grid, see equation (7) in Philippe's original paper. For the U-flux, the velocity conversion gives the deg2m * math.cos(rad * y) factor and the edge gives the deg2m factor, for the V-flux, the velocity conversion gives the deg2m factor and the edge gives the deg2m * math.cos(rad * y). Thus the total conversion factor for both the U and V flux is meshJac = deg2m * deg2m * math.cos(rad * y)

erikvansebille · 2026-01-09T12:55:25Z

Thanks for this explanation, @MeikeBos. Makes sense, yes!

It's a relief that we now also understand why the cosine-factor is also needed for the meridional velocity 😃

Following the discussion around zonal and meridional spherical conversion in #2455

Reverting cgrid velocity conversion back to v3 implementation

46ffe34

github-project-automation bot added this to Parcels development Jan 8, 2026

github-project-automation bot moved this to Backlog in Parcels development Jan 8, 2026

Adding unit test for meridional advection on curvilinear grid

5436136

To test that the spherical mesh conversion leads to the expected meridional velocity

erikvansebille marked this pull request as ready for review January 9, 2026 12:55

erikvansebille requested a review from VeckoTheGecko January 9, 2026 12:55

erikvansebille mentioned this pull request Jan 9, 2026

Removing UnitConverters from v4 #2459

Merged

erikvansebille added a commit that referenced this pull request Jan 12, 2026

Making unit tests of zonal and meridional flow stricter

49a0a82

Following the discussion around zonal and meridional spherical conversion in #2455

erikvansebille mentioned this pull request Jan 12, 2026

Implementing convert.nemo_to_sgrid() #2460

Merged

VeckoTheGecko approved these changes Jan 12, 2026

View reviewed changes

github-project-automation bot moved this from Backlog to Ready in Parcels development Jan 12, 2026

erikvansebille merged commit d5009e1 into fieldset_from_nemo Jan 12, 2026
11 checks passed

erikvansebille deleted the unitconversion_for_cgridvelocity branch January 12, 2026 13:17

github-project-automation bot moved this from Ready to Done in Parcels development Jan 12, 2026

erikvansebille mentioned this pull request Jan 13, 2026

Adding XLinear_Velocity interpolator #2461

Merged

3 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Fixing unit conversion for C-grid velocity interpolation #2455

Fixing unit conversion for C-grid velocity interpolation #2455

Uh oh!

erikvansebille commented Jan 8, 2026 •

edited

Loading

Uh oh!

erikvansebille commented Jan 9, 2026

Uh oh!

MeikeBos commented Jan 9, 2026

Uh oh!

erikvansebille commented Jan 9, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

	rad = np.pi / 180.0
	deg2m = 1852 * 60.0
	if applyConversion:
	meshJac = (deg2m * deg2m * math.cos(rad * y)) if grid.mesh == "spherical" else 1
	else:
	meshJac = deg2m if grid.mesh == "spherical" else 1

	jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) * meshJac

	u = (
	(-(1 - eta) * U - (1 - xsi) * V) * px[0]
	+ ((1 - eta) * U - xsi * V) * px[1]
	+ (eta * U + xsi * V) * px[2]
	+ (-eta * U + (1 - xsi) * V) * px[3]
	) / jac
	v = (
	(-(1 - eta) * U - (1 - xsi) * V) * py[0]
	+ ((1 - eta) * U - xsi * V) * py[1]
	+ (eta * U + xsi * V) * py[2]
	+ (-eta * U + (1 - xsi) * V) * py[3]
	) / jac

Fixing unit conversion for C-grid velocity interpolation #2455

Fixing unit conversion for C-grid velocity interpolation #2455

Uh oh!

Conversation

erikvansebille commented Jan 8, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

erikvansebille commented Jan 9, 2026

Uh oh!

MeikeBos commented Jan 9, 2026

Uh oh!

erikvansebille commented Jan 9, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

erikvansebille commented Jan 8, 2026 •

edited

Loading